"""
CRTM Fortran Wrapper Bindings for Julia

Provides Julia ccall interface to the CRTM Fortran wrapper library (libcrtm_wrapper.so).
Enables zero-copy satellite radiance forward and adjoint operators.

# Functions
- `crtm_init_julia` - Initialize CRTM with sensor configuration
- `crtm_forward_julia` - Forward radiative transfer operator H(x)
- `crtm_adjoint_julia` - Adjoint operator H^T(y)
- `crtm_destroy_julia` - Cleanup CRTM resources
- `crtm_version_julia` - Get CRTM version string

# Example
```julia
using GSICoreAnalysis.FortranInterface

# Initialize CRTM
status = crtm_init_julia("amsua_n15", "/path/to/coefficients", [1,2,3,4,5])

# Forward operator
temp = rand(Float32, 32, 10) .+ 250.0f0
pres = rand(Float32, 32, 10) .* 1000.0f0 .+ 100.0f0
humid = rand(Float32, 32, 10) .* 0.01f0
tskin = rand(Float32, 10) .+ 280.0f0
wind = rand(Float32, 10) .* 5.0f0
zenith = rand(Float32, 10) .* 60.0f0
azimuth = zeros(Float32, 10)

tb = crtm_forward_julia(temp, pres, humid, tskin, wind, zenith, azimuth, [1,2,3,4,5])

# Adjoint operator
innovation = rand(Float32, 5, 10) .* 2.0f0
temp_ad, humid_ad = crtm_adjoint_julia(innovation, temp, pres, humid,
                                       tskin, wind, zenith, azimuth, [1,2,3,4,5])
```
"""

# C-compatible structure definitions matching Fortran ISO_C_BINDING types
struct CAtmosphere
    n_layers::Int32
    n_profiles::Int32
    temperature::Ptr{Float32}
    pressure::Ptr{Float32}
    humidity::Ptr{Float32}
end

struct CSurface
    n_profiles::Int32
    skin_temperature::Ptr{Float32}
    wind_speed::Ptr{Float32}
    surface_type::Int32
end

struct CGeometry
    n_profiles::Int32
    zenith_angle::Ptr{Float32}
    azimuth_angle::Ptr{Float32}
end

# Path to compiled CRTM wrapper library
const LIBCRTM_WRAPPER = joinpath(@__DIR__, "../../deps/libcrtm_wrapper.so")

"""
    crtm_init_julia(sensor_id::String, coef_path::String, channels::Vector{Int32}) -> Int32

Initialize CRTM with sensor configuration.

# Arguments
- `sensor_id`: Sensor identifier (e.g., "amsua_n15", "iasi_metop-a")
- `coef_path`: Path to CRTM coefficient files
- `channels`: Channel numbers to use (1-indexed)

# Returns
- Status code (0 = success, non-zero = error)
"""
function crtm_init_julia(sensor_id::String, coef_path::String, channels::Vector{Int32})
    sensor_bytes = Vector{UInt8}(sensor_id)
    coef_bytes = Vector{UInt8}(coef_path)
    n_channels = Int32(length(channels))
    status = Ref{Int32}(0)

    ccall((:crtm_init_c, LIBCRTM_WRAPPER),
          Cvoid,
          (Ptr{UInt8}, Int32, Ptr{UInt8}, Int32, Int32, Ptr{Int32}, Ref{Int32}),
          sensor_bytes, Int32(length(sensor_bytes)),
          coef_bytes, Int32(length(coef_bytes)),
          n_channels, channels, status)

    return status[]
end

"""
    crtm_forward_julia(temperature, pressure, humidity, skin_temp, wind_speed,
                       zenith_angle, azimuth_angle, channels) -> Matrix{Float32}

CRTM forward operator: atmospheric state → brightness temperatures.

# Arguments
- `temperature`: Temperature profile [K] (n_layers × n_profiles)
- `pressure`: Pressure profile [hPa] (n_layers × n_profiles)
- `humidity`: Specific humidity [kg/kg] (n_layers × n_profiles)
- `skin_temp`: Surface skin temperature [K] (n_profiles)
- `wind_speed`: Surface wind speed [m/s] (n_profiles)
- `zenith_angle`: Satellite zenith angle [degrees] (n_profiles)
- `azimuth_angle`: Satellite azimuth angle [degrees] (n_profiles)
- `channels`: Channel numbers to simulate (1-indexed)

# Returns
- `brightness_temp`: Simulated brightness temperatures [K] (n_channels × n_profiles)

# Notes
- All arrays must be Float32
- Column-major layout (Julia/Fortran compatible)
- Zero-copy: arrays passed by pointer
"""
function crtm_forward_julia(temperature::Matrix{Float32},
                            pressure::Matrix{Float32},
                            humidity::Matrix{Float32},
                            skin_temp::Vector{Float32},
                            wind_speed::Vector{Float32},
                            zenith_angle::Vector{Float32},
                            azimuth_angle::Vector{Float32},
                            channels::Vector{Int32})

    # Validate dimensions
    nlev, nprof = size(temperature)
    @assert size(pressure) == (nlev, nprof) "Pressure size mismatch"
    @assert size(humidity) == (nlev, nprof) "Humidity size mismatch"
    @assert length(skin_temp) == nprof "Surface temperature size mismatch"
    @assert length(wind_speed) == nprof "Wind speed size mismatch"
    @assert length(zenith_angle) == nprof "Zenith angle size mismatch"
    @assert length(azimuth_angle) == nprof "Azimuth angle size mismatch"

    nchans = Int32(length(channels))

    # Create C-compatible structures
    atm = CAtmosphere(
        Int32(nlev),
        Int32(nprof),
        pointer(temperature),
        pointer(pressure),
        pointer(humidity)
    )

    sfc = CSurface(
        Int32(nprof),
        pointer(skin_temp),
        pointer(wind_speed),
        Int32(0)  # surface_type (0 = ocean default)
    )

    geom = CGeometry(
        Int32(nprof),
        pointer(zenith_angle),
        pointer(azimuth_angle)
    )

    # Allocate output
    brightness_temp = zeros(Float32, nchans, nprof)
    status = Ref{Int32}(0)

    # Call Fortran wrapper
    ccall((:crtm_forward_c, LIBCRTM_WRAPPER),
          Cvoid,
          (Ref{CAtmosphere}, Ref{CSurface}, Ref{CGeometry},
           Int32, Ptr{Int32}, Ptr{Float32}, Ref{Int32}),
          Ref(atm), Ref(sfc), Ref(geom),
          nchans, channels, brightness_temp, status)

    if status[] != 0
        error("CRTM forward operator failed with status $(status[])")
    end

    return brightness_temp
end

"""
    crtm_adjoint_julia(tb_innovation, temperature, pressure, humidity,
                       skin_temp, wind_speed, zenith_angle, azimuth_angle,
                       channels) -> (temp_adjoint, humid_adjoint)

CRTM adjoint operator: brightness temperature innovations → state-space gradients.

# Arguments
- `tb_innovation`: Brightness temperature innovations [K] (n_channels × n_profiles)
- `temperature`: Temperature profile [K] (n_layers × n_profiles)
- `pressure`: Pressure profile [hPa] (n_layers × n_profiles)
- `humidity`: Specific humidity [kg/kg] (n_layers × n_profiles)
- `skin_temp`: Surface skin temperature [K] (n_profiles)
- `wind_speed`: Surface wind speed [m/s] (n_profiles)
- `zenith_angle`: Satellite zenith angle [degrees] (n_profiles)
- `azimuth_angle`: Satellite azimuth angle [degrees] (n_profiles)
- `channels`: Channel numbers (1-indexed)

# Returns
- `temp_adjoint`: Adjoint forcing for temperature [K] (n_layers × n_profiles)
- `humid_adjoint`: Adjoint forcing for humidity [kg/kg] (n_layers × n_profiles)

# Notes
- Adjoint consistency: <H(x), y> = <x, H^T(y)>
- Used in variational minimization gradient computation
"""
function crtm_adjoint_julia(tb_innovation::Matrix{Float32},
                            temperature::Matrix{Float32},
                            pressure::Matrix{Float32},
                            humidity::Matrix{Float32},
                            skin_temp::Vector{Float32},
                            wind_speed::Vector{Float32},
                            zenith_angle::Vector{Float32},
                            azimuth_angle::Vector{Float32},
                            channels::Vector{Int32})

    # Validate dimensions
    nchans, nprof = size(tb_innovation)
    nlev, nprof2 = size(temperature)
    @assert nprof == nprof2 "Profile count mismatch"
    @assert Int32(length(channels)) == nchans "Channel count mismatch"
    @assert size(pressure) == (nlev, nprof) "Pressure size mismatch"
    @assert size(humidity) == (nlev, nprof) "Humidity size mismatch"

    # Create C-compatible structures
    atm = CAtmosphere(
        Int32(nlev),
        Int32(nprof),
        pointer(temperature),
        pointer(pressure),
        pointer(humidity)
    )

    sfc = CSurface(
        Int32(nprof),
        pointer(skin_temp),
        pointer(wind_speed),
        Int32(0)
    )

    geom = CGeometry(
        Int32(nprof),
        pointer(zenith_angle),
        pointer(azimuth_angle)
    )

    # Allocate output
    temp_ad = zeros(Float32, nlev, nprof)
    humid_ad = zeros(Float32, nlev, nprof)
    status = Ref{Int32}(0)

    # Call Fortran wrapper
    ccall((:crtm_adjoint_c, LIBCRTM_WRAPPER),
          Cvoid,
          (Ptr{Float32}, Int32, Int32, Int32,
           Ref{CAtmosphere}, Ref{CSurface}, Ref{CGeometry},
           Ptr{Int32}, Ptr{Float32}, Ptr{Float32}, Ref{Int32}),
          tb_innovation, Int32(nchans), Int32(nprof), Int32(nlev),
          Ref(atm), Ref(sfc), Ref(geom),
          channels, temp_ad, humid_ad, status)

    if status[] != 0
        error("CRTM adjoint operator failed with status $(status[])")
    end

    return temp_ad, humid_ad
end

"""
    crtm_version_julia() -> String

Get CRTM wrapper version string.
"""
function crtm_version_julia()
    version_len = 64
    version_str = zeros(UInt8, version_len)

    ccall((:crtm_version_c, LIBCRTM_WRAPPER),
          Cvoid,
          (Int32, Ptr{UInt8}),
          Int32(version_len), version_str)

    # Find first zero byte, or use full string if none found
    zero_idx = findfirst(iszero, version_str)
    if isnothing(zero_idx)
        return String(version_str)
    else
        return String(version_str[1:zero_idx-1])
    end
end

"""
    crtm_destroy_julia() -> Int32

Cleanup CRTM resources.

# Returns
- Status code (0 = success)
"""
function crtm_destroy_julia()
    status = Ref{Int32}(0)

    ccall((:crtm_destroy_c, LIBCRTM_WRAPPER),
          Cvoid,
          (Ref{Int32},),
          status)

    return status[]
end

# Export functions
export crtm_init_julia, crtm_forward_julia, crtm_adjoint_julia
export crtm_version_julia, crtm_destroy_julia
export CAtmosphere, CSurface, CGeometry
